home *** CD-ROM | disk | FTP | other *** search
/ FishMarket 1.0 / FishMarket v1.0.iso / fishies / 376-400 / disk_386 / xlispstat / src2.lzh / XLisp-Stat / derivatives.c < prev    next >
C/C++ Source or Header  |  1990-10-02  |  2KB  |  79 lines
  1. /* derivatives - for Bayes code in XLISP-STAT and S                    */
  2. /* XLISP-STAT 2.1 Copyright (c) 1990, by Luke Tierney                  */
  3. /* Additions to Xlisp 2.1, Copyright (c) 1989 by David Michael Betz    */
  4. /* You may give out copies of this software; for conditions see the    */
  5. /* file COPYING included with this distribution.                       */
  6.  
  7. #include "xlisp.h"
  8. #include "osdef.h"
  9. #ifdef ANSI
  10. #include "xlsproto.h"
  11. #else
  12. #include "xlsfun.h"
  13. #endif ANSI
  14.  
  15. # define nil 0L
  16. # define FALSE 0
  17. # define TRUE 1
  18.  
  19. void numergrad(n, x, grad, fsum, ffun, h, typx)
  20.      int n;
  21.      RVector x, grad, fsum, typx;
  22. #ifdef ANSI
  23.      int (*ffun)(RVector,double *,int,int);
  24. #else
  25.      int (*ffun)();
  26. #endif ANSI
  27.      double h;
  28. {
  29.   int i;
  30.   double old_xi, f1, f2, hi;
  31.  
  32.   for (i = 0; i < n; i++) {
  33.     old_xi = x[i];
  34.     hi = (typx != nil) ? typx[i] * h : h;
  35.     x[i] = old_xi + hi;
  36.     (*ffun)(x, &f1, nil, nil);
  37.     x[i] = old_xi - hi;
  38.     (*ffun)(x, &f2, nil, nil);
  39.     x[i] = old_xi;
  40.     grad[i] = (f1 - f2) / (2.0 * hi);
  41.     fsum[i] = f1 + f2;
  42.   }
  43. }
  44.  
  45. void numerhess(n, x, hess, f, fsum, ffun, h, typx)
  46.      int n;
  47.      RVector x, fsum, typx;
  48.      RMatrix hess;
  49. #ifdef ANSI
  50.      int (*ffun)(RVector,double *,int,int);
  51. #else
  52.      int (*ffun)();
  53. #endif ANSI
  54.      double h, f;
  55. {
  56.   int i, j;
  57.   double old_xi, old_xj, f1, f2, hi, hj;
  58.  
  59.   for (i = 0; i < n; i++) {
  60.     hi = (typx != nil) ? typx[i] * h : h;
  61.     hess[i][i] = (fsum[i] - 2 * f) / (hi * hi);
  62.     for (j = i + 1; j < n; j++) {
  63.       hj = (typx != nil) ? typx[j] * h : h;
  64.       old_xi = x[i];
  65.       old_xj = x[j];
  66.       x[i] = old_xi + hi;
  67.       x[j] = old_xj + hj;
  68.       (*ffun)(x, &f1, nil, nil);
  69.       x[i] = old_xi - hi;
  70.       x[j] = old_xj - hj;
  71.       (*ffun)(x, &f2, nil, nil);
  72.       x[i] = old_xi;
  73.       x[j] = old_xj;
  74.       hess[i][j] = (2 * f + f1 + f2 - fsum[i] - fsum[j]) / (2.0 * hi * hj);
  75.       hess[j][i] = hess[i][j];
  76.     }
  77.   }
  78. }
  79.